function [derivs] = grad_3_44(hyps, func, n, n_class, X, y, approxF)
% a function evaluating equation 3.44 and derivatives, of R&W
% by Mark Norrish, 2011
% note: actually uses a bunch of other equations, esp. for derivations;
% section 5.5.1 in particular

% hyps: as a vector; [ h_class1 ; ... ; h_class_c ]

dim = length(hyps) / n_class;
Hyps = reshape(hyps, dim, n_class);

bigK = zeros(n*n_class); K = zeros(n,n,n_class); sigma_noise = 1e-7;
% NTS: rewrite without bigK, it's just embarrassing itself as is
for c = 1:n_class
  K(:,:,c) = func(Hyps(:,c), X, X) + sigma_noise*eye(n);
  bigK(1+(c-1)*n:c*n,1+(c-1)*n:c*n) = K(:,:,c);
end

if nargin <= 6
  f = alg_3_3(n, n_class, K, y);
else
  f = alg_3_3(n, n_class, K, y, approxF);
end
F = reshape(f, n, n_class); % it's SO much easier to think with
expsum = repmat(sum(exp(F)')',n_class,1);
pi = exp(f)./expsum;
bigPi = zeros(n_class*n,n);
for i = 1:n_class
  bigPi((i-1)*n+1:i*n,:) = diag(pi((i-1)*n+1:i*n));
end
W = diag(pi) - bigPi * bigPi';
Pi = reshape(pi,n,n_class);
%[lW, uW] = lu(W);
derivs = zeros(size(hyps));

expsum = expsum(1:n);
dpyfdf = reshape(y, n, n_class) - exp(F) ./ repmat(expsum, 1, n_class); % d log p(y|f) / d f
% NTS: rewrite to need fewer reshapes
IKW = (1 + sigma_noise) * eye(n*n_class) + bigK * W;
IKWK = (IKW \ bigK); % for eqn 5.23; NTS: this can be computed O(c^2)-faster if done blockwise
WIKW = W / IKW;
for c = 1:n_class
  crn = 1+(c-1)*n:c*n; % class-range: i.e. the class-appropriate block
  alpha = K(:,:,c) \ F(:,c);
  thisIKW = IKW(crn,crn); % consider storing the inverse of this: faster if less stable
  equalities = repmat(c == 1:n_class,n_class,1);
  del_W = zeros(n_class * n);
  for i = 1:n
    del_W = 2 * Pi(i,c) * Pi(i,:)' * Pi(i,:) - (equalities + equalities') .* (Pi(i,:)'*Pi(i,:)) - diag(Pi(i,:)*Pi(i,c));
    del_W(c,c) += Pi(i,c);
    %for c1 = 1:n
    %  del_W((c-1)*n+i,(c1-1)*n + i) = 
    %end
    z = zeros(n); z(i,i) = 1;
    del_W = kron(del_W, z);
    del_logq(i) = trace(IKWK * del_W);
  end
  for hyp_id = 1:dim
    h = (c-1)*dim+hyp_id;
    del_K = func(Hyps(:, c), X, X, hyp_id);
    del_f = thisIKW \ del_K * dpyfdf(:, c); % d fhat / d theta_i; by eq 5.24
    del_logq = zeros(n,1); % dlog(q)/d(fhat_i^c)
    
    % eq 5.22:
    %derivs(h) = 0.5 * alpha' * del_K * alpha;% - 0.5 * trace(W(crn,crn) / thisIKW * del_K) - 0.5 * sum(del_logq .* del_f);
    derivs(h) = -0.5 * trace(WIKW(crn,crn) * del_K) - 0.5 * sum(del_logq .* del_f);
  end
end

derivs = -derivs; % because it's negloglik not loglik